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Abstract 

We propose a new spectral Lagrangian based deterministic solver for the non-linear 
Boltzmann Transport Equation for Variable Hard Potential (VHP) collision kernels 
with conservative or non-conservative binary interactions. The method is based on 
symmetries of the Fourier transform of the collision integral, where the complex- 
ity in its computing is reduced to a separate integral over the unit sphere 5^. In 
addition, the conservation of moments is enforced by Lagrangian constraints. The 
resulting scheme, implemented in free space is very versatile and adjusts in a very 
simple manner, to several cases that involve energy dissipation due to local micro- 
reversibility (inelastic interactions) or elastic model of slowing down process. Our 
simulations are benchmarked with the available exact self-similar solutions, exact 
moment equations and analytical estimates for homogeneous Boltzmann equation 
for both elastic and inelastic VHP interactions. Benchmarking of the simulations 
involves the selection of a time self-similar rescaling of the numerical distribution 
function which is performed using the continuous spectrum of the equation for 
Maxwell molecules as studied first in [13] and generalized to a wide range of related 
models in [12]. The method also produces accurate results in the case of inelastic 
diffusive Boltzmann equations for hard-spheres (inelastic collisions under thermal 
bath), where overpopulated non-Gaussian exponential tails have been conjectured 
in computations by stochastic methods in [49; 26; 46; 35] and rigourously proven in 
[34] and [15]. 
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1 Introduction 



In a microscopic description of a rarefied gas, all particles are assumed to be 
traveling in a straight line with a fixed velocity until they enter into a colli- 
sion. In such dilute flows, binary collisions are often assumed to be the main 
mechanism of particle interactions. The statistical effect of such collisions can 
be modeled by collision terms of the Boltzmann or Enskog transport equa- 
tion type, where the kinetic dynamics of the gas are subject to the molecular 
chaos assumption. The nature of these interactions could be elastic, inelas- 
tic or coalescing. They could either be isotropic or anisotropic, depending on 
their collision rates as a function of the scattering angle. In addition, collisions 
are described in terms of inter-particle potentials and the rate of collisions is 
usually modeled as product of power laws for the relative speed and the differ- 
ential cross section, at the time of the interaction. When the rate of collisions 
is independent of the relative speed, the interaction is referred to as of Maxwell 
type. When it corresponds to relative speed to a positive power less than unity, 
they are referred to as Variable Hard Potentials (VHP) and when the rate of 
collisions is proportional to the relative speed, it is referred to as hard spheres. 

The Boltzmann Transport Equation (an integro-differential transport equa- 
tion) describes the evolution of a single point probability distribution function 
f{x,v,t) which is defined as the probability of finding a particle at position 
X with velocity (kinetic) v at time t. The mathematical and computational 
difficulties associated to the Boltzmann equation are due to the non local - 
non linear nature of the collision operator, which is usually modeled as a multi 
linear integral form in c?-dimensional velocity space and unit sphere S'^~^. 



From the computational point of view, of the well-known and well-studied 
methods developed in order to solve this equation is an stochastic based 
method called "Direct Simulation Monte-Carlo" (DSMC) developed initially 
by Bird [5] and Nanbu 48 1 and more recently by js^; 55]. This method is 
usually employed as an alternative to hydrodynamic solvers to model the evo- 
lution of moments or hydrodynamic quantities. In particular, this method have 
been shown to converge to the solution of the classical Boltzmann equation 
in the case of mono atomic rarefied gases [57]. One of the main drawbacks of 
such methods is the inherent statistical fluctuations in the numerical results, 
which becomes very expensive or unreliable in presence of non-stationary flows 
or non equilibrium statistical states, where more information is desired about 
the evolving probability distribution. Currently, there is extensive work from 
Rjasanow and Wagner jH] and references therein, to determine accurately 
the high- velocity tail behavior of the distribution functions from DSMC data. 
Implementations for micro irreversible interactions such as inelastic collisions 



have been carefully studied in [35 
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In contrast, a deterministic method computes approximations of the prob- 
abihty distribution function using the Boltzmann equation, as well as ap- 
proximations to the observables like density, momentum, energy, etc.,. There 
are currently two deterministic approaches to the computations of non-linear 
Boltzmann, one is the well known discrete velocity models and the second 
a spectral based method, both implemented for simulations of elastic inter- 
actions i.e. energy conservative evolution. Discrete velocity models were de- 
veloped by Broadwell [2(J] and mathematically studied by Illner, Cabannes, 
Kawashima among many authors 4l|; |42; |21|. More recently these models 



have been studied for many other applications on kinetic elastic theory in 
0;[24; 44; 59; 39 1. These models have not adapted to inelastic coUisional prob- 
lems up to this point according to our best knowledge. 



Spectral based models, which are the ones of our choice in this work, have 



been developed by Pareschi, Gabetta and Toscani [32|, and later by Bobylev 



and Rjasanow [l7j and Pareschi and Russo l52|. These methods are supported 
by the ground breaking work of Bobylev |4|] using the Fourier Transformed 
Boltzmann Equation to analyze its solutions in the case of Maxwell type 
of interactions. After the introduction of the inelastic Boltzmann equation 
for Maxwell type interactions and the use of the Fourier transform for its 
analysis by Bobylev, Carrillo and one of the authors here [gI, the spectral 
based approach is becoming the most suitable tool to deal with determin- 
istic computations of kinetic models associated with Boltzmann non-linear 
binary coUisional integral, both for elastic or inelastic interactions. More re- 
cent implementations of spectral methods for the non-linear Boltzmann are 
due to Bobylev and Rjasanow 
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who developed a method using the Fast 
Fourier Transform (FFT) for Maxwell type of interactions and then for Hard- 
Sphere interactions [isj using generalized Radon and X-ray transforms via 



FFT. Simultaneously, L. Pareschi and B. Perthame [5l[ developed similar 
scheme using FFT for Maxwell type of interactions. Later, I. Ibragimov and 
S. Rjasanow [40| developed a numerical method to solve the space homoge- 
neous Boltzmann Equation on a uniform grid for a Variable Hard Potential 
interactions with elastic collisions. This particular work has been a great in- 
spiration for the current work and was one of the first initiating steps in the 
direction of a new numerical method. 



We mention that, most recently, Filbet and Russo [27[, [28| implemented a 
method to solve the space inhomogeneous Boltzmann equation using the pre- 
viously developed spectral methods in 52|;|51|- Afore mentioned work in devel- 
oping deterministic solvers for non-linear BTE have been restricted to elastic, 
conservative interactions. Finally, Mouhout and Pareschi 47| are currently 
studying the approximation properties of the schemes. Part of the difficulties 
in their strategy arises from the constraint that the numerical solution has to 
satisfy conservation of the initial mass. To this end, the authors propose the 
use of a periodic representation of the distribution function to avoid aliasing. 
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There is no conservation of momentum and energy in 28|] , 27|] and 47| . Both 



methods ( [28|, [27[, [47j), which are developed in 2 and 3 dimensions, do not 
guarantee the positivity of the solution due to the fact that the truncation of 
the velocity domain combined with the Fourier method makes the distribution 
function negative at times. This last shortcoming of the spectral approach re- 
mains in our proposed technique; however we are able to handle conservation 
in a very natural way by means of Lagrange multipliers. We also want to credit 
an unpublished calculation of V. Panferov and S. Rjasanow ^] who wrote 
a method to calculate the particle distribution function for inelastic collisions 
in the case of hard spheres, but there were no numerical results to corroborate 
the efficiency of the method. Our proposed approach is slightly different and 
it takes a less number of operations to compute the collision integral. 

Our current approach, based on a modified version of the work in and 
4o| . works for elastic or inelastic collisions and energy dissipative non-linear 
Boltzmann type models for variable hard potentials. We do not use periodic 
representations for the distribution function. The only restriction of the cur- 
rent method is that it requires that the distribution function at any time 
step be Fourier transformable. The required conservation properties of the 
distribution function are enforced through a Lagrange multiplier constrained 
optimization problem with the desired conservation quantities set as the con- 
straints. Such corrections to the distribution function to make it conservative 
are very small but crucial for the evolution of the probability distribution 
function according to the Boltzmann equation. 

This Lagrange optimization problem gives the freedom of not conserving the 
energy, independent of the collision mechanism, as long momentum is con- 
served. Such a technique plays a major role as it gives the option of com- 
puting energy dissipative solutions by just eliminating one constraint in the 
corresponding optimization problem. The current method can be easily im- 
plemented in any dimension. A novel aspect of the presented approach here 
lays on a new method that uses the Fourier Transform as a tool to simplify 
the computation of the collision operator that works, both for elastic and 
inelastic collisions. It is based on an integral representation of the Fourier 



Transform of the collision kernel as used in [17[. If is the number of dis- 
cretizations in one direction of the velocity domain in (i- dimensions, the total 
number of operations required to solve for the collision integral are of the 
order of N'^Hog{N) + 0{N'^'^). And this number of operations remains the 
same for elastic/ inelastic, isotropic/ anisotropic VHP type of interactions. 
However, when the differential cross section is independent of the scattering 
angle, the integral representation kernel is further reduced by an exact closed 
integrated form that is used to save in computational number of operations 
to 0{N'^log{N)). This reduction is possible when computing hard spheres in 
d+2 dimensions or Maxwell type models in 2-dimensions. Nevertheless, the 
method can be employed without much changes for the other case. In partic- 
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ular the method becomes 0{P'^~^ N'^log{N)), where P, the number of each 
angular discretizations is expected to be much smaller than used for energy 
discretizations. Such reduction in number of operations was also reported in 



28l | with 0{Nlog{N)) number of operations, where the authors are assuming 



N to be the total number of discretizations in the ci-dimensional space (i.e. 
our A^'^ and P of order of unity). 



Our numerical study is performed for several examples of well establish behav- 
ior associated to solutions of energy dissipative space homogeneous collisional 
models under heating sources that secure existence of stationary states with 
positive and finite energy. We shall consider heating sources corresponding 
to randomly heated inelastic particles in a heat bath, with and without fric- 
tion; elastic or inelastic collisional forms with anti-divergence terms due to 
dynamically (self-similar) energy scaled solutions 3J; [iSj and a particularly 
interesting example of inelastic collisions added to a slow down linear process 
that can be derived as a weakly coupled heavy/light binary mixture. On this 
particular case, when Maxwell type interactions are considered, it is shown 



that [13|; [IJ; [12], on one hand dynamically energy scaled solutions exist, they 



have a close, explicit formula in Fourier space for a particular choice of pa- 
rameters and their corresponding anti Fourier transform in probability space 
exhibit a singularity at the origin and power law high energy tails, while re- 
maining integrable and with finite energy. On the other hand they are stable 
within a large class of initial states. We used this particular example to bench- 
mark our computations by spectral methods by comparing the dynamically 
scaled computed solutions to the explicit one self similar one. 



Convergence and error results of the Fourier Transform Lagrangian method, 
locally in time, are currently being developed [36j, and it is expected that the 



proposed spectral approximation of the free space problem will have optimal 
algorithm complexity using the non-equispaced FFT as obtained by Greengard 
and Lin [s^ for spectral approximation of the free space heat kernel. 



Implementation of the space inhomogeneous case are also currently being con- 
sidered. The spectral-Lagrangian scheme methodology proposed here can be 
extended to cases of Pareto tails, opinion dynamics and N player games, where 
the evolution and asymptotic behavior of probabilities are studied in Fourier 



space as well. |53|; [12 



The paper is organized as follows. In section 2, some preliminaries and de- 
scription of the various approximated models associated with the elastic or 
inelastic Boltzmann equation are presented. In section 3, the actual numerical 
method is discussed with a small discussion on its discretization. In section 
4, the special case of spatially homogeneous collisional model for a slow down 
process derived from a weakly coupled binary problem with isotropic elastic 
Maxwell type interactions is considered wherein an explicit solution is derived 
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and shown to have power-hke tails in some particular cases corresponding to 
a cold thermostat problem. Section 5 deals with the numerical results and 
examples. Finally in section 6, direction of future work is proposed along with 
a summary of the proposed numerical method. 



2 Preliminaries 



The initial value problem associated to space homogeneous Boltzmann Trans- 
port Equation modeling the statistical (kinetic) evolution of a single point 
probability distribution function f{v,t) for Variable Hard Potential (VHP) 
interactions is given by 



^/(^,t) = Qif,f)iv,t) 



Wi'v, t)f{'w, t) - f{v, t)f{w, t)] B{\u\,fi) dadw 

f{v,0) = fo{v), (2.1) 

where the initial probability distribution fo{v) is assumed integrable and 
Jfs = is Jacobian of post with respect to pre coUisional velocities which 

depend the local energy dissipation The problem may or may not have fi- 
nite initial energy Sq = /^^ fo{v)\v\'^dv and the velocity interaction law, written 
in center of mass and relative velocity coordinates is 



(2.2) 



u = V — w : the relative velocity 

v' = V + ^{\u\a — u), w' = w — ^{\u\a — u) , 
11 . (J 

a = cos(6) = —. — j- : the cosine of the scattering angle , 

\u\ 

B{\u\,iJ,) = \u\^b{cose) withO<A<l, 
/ b(cos 6) sin'^^'^ 9 dO < K : Grad cut-off assumption 
1 + e 

(3 = : the energy dissipation parameter , 

2 

where the parameter e G [0, 1] is the restitution coefficient corresponding from 
sticky to elastic interactions, where Jg = Ji = 1. 

We denote by 'v and 'w the pre-collision velocities corresponding to v and w. 
In the case of micro-reversible (elastic) collisions one can replace 'v and 'w with 
v' and w' respectively in the integral part of (2.1). We assume the differential 
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cross section function b{j^-) is integrable with respect to the post-colhsional 
specular reflection direction a in the d—1 dimensional sphere, referred as the 
Grad cut-off assumption, and that 6(cos 6) is renormalized such that 



where the constant uJd-2 is the measure of the d — 2 dimensional sphere and 
the corresponding scattering angle is 6 is defined by cos 6* = 

The parameter A regulates the collision frequency as a function of the rel- 
ative speed \u\. It accounts for inter particle potentials defining the coUisional 
kernel and they are referred to as Variable Hard Potentials (VHP) whenever 
< A < 1, Maxwell Molecules type interactions (MM) for A = and Hard 
Spheres (HS) for A = 1. The Variable Hard Potential collision kernel then 
takes the following general form: 



with Cx{a) = ■^b{9),X = for Maxwell type of interactions; C\{(r) = 

2 

^,A = 1 for Hard Spheres. In addition, if C\{a) is independent of the 
scattering angle we call the interactions isotropic. Otherwise we referred to 
them as anisotropic Variable Hard Potential interactions. 

For classical case of elastic collisions, it has been established that the Cauchy 
problem for the space homogeneous Boltzmann equation has a unique solu- 
tion in the class of integrable functions with finite energy (i.e. C^i^LliW^))), 
it is regular if initially so, and /(.,t) converges in LKW^) to the Maxwellian 
distribution Mpy^g{v) associated to the d + 2- moments of the initial state 
f{v,0) = fo{v) e L\{M'^). In addition, if the initial state has Maxwellian de- 
cay, this property will remain with a Maxwellian decay globally bounded in 
time ([s^]), as well as all derivatives if initial so (see [H). 

Depending on their nature, collisions either conserve density, momentum and 
energy (elastic) or density and momentum (inelastic) or density (elastic - lin- 
ear Boltzmann operator), depending on the number of collision invariants the 
operator Q{f, f){t, v) has. In the case of the classical Boltzmann equation for 
rarefied (elastic) mono-atomic gases, the collision invariants are exactly d + 2, 
that is, according to the Boltzmann theorem, the number of polynomials in 
velocity space v that generate 0(f) = A + B ■ v + C|v|^, with C < 0. In 
particular, one obtains the following conserved quantities 




(2.3) 



B{\u\,fi) = Cx{a)\u 



X 



(2.4) 
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density p{t) = / f{v^t)dv 

momentum m{t) = / vf{v,t)dv (2-5) 



vet. 

1 ' . ,2. 



energy £{t) = — -— / |f| f{v,t)dv. 



Of significant interest from the statistical view point are the evolution of mo- 
ments or observables, at all orders. They are defined by the djTiamics of the 
corresponding time evolution equation for the velocity averages, given by 



^M,(t) = / f{v,t)v®^dv = I QifJ)iv,t)v®=dv , (2.6) 



where, v®^ = the standard symmetric tensor product of v with itself, j times. 
Thus, according to (2.6), for the classical elastic Boltzmann equation, the first 
d + 2 moments are conserved, meaning, Mj{t) = Mqj = / ^ fQ{v)v®^dv for 

j = 0, 1; and S{t) = tr(M2)(t) = Sq = J fo(v)\v\'^dv. Higher order moments 
or observables of interest are 



Momentum Flow M2{t) = I vv^f{v,t)dv 

Energy Flow r(t) = — -— / vlvl"^ f(v,t)dv 

2p(t) 4 I I M , ; 



Bulk Velocity V{t) = ^ (2.7) 
Internal Energy £{t) = —{tr{M2)-p\V\'^) 



Temperature T(t) 



with k— Boltzmann constant. 



2p 
28 jt) 

kd 



We finally point out that, in the case of Maxwell molecules (A = 0), it is 
possible to write recursion formulas for higher order moments of all orders (j^ 
for the elastic case, and [6(] in the inelastic case) which, in the particular case 
of isotropic solutions depending only on |f P/2, take the form 
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m„(t) = / \v\''' f{v,t)dv = e-^"*m„(0) + 



n-l 



with (2.5 
1 

+ 



A„ = i-^[/3'"+E(i-/5r] 



Baik.n-k) = p^'' C s^il- (3(2- (3)sY-^ds, 
Jo 

for 72 > 1, < /? < 1, where Aq = 0, mo(t) = 1, and m„,(0) = / |f fo{v)dv. 



2. 1 Boltzmann collisional models with heating sources 



A colhsional model associated to the space homogeneous Boltzmann transport 
equation (2.1) with grad cutoff assumption (2.2), can be modified in order to 
accommodate for an energy or 'heat source' like term Q{f{t,v)), where Q 
is a differential or integral operator. In these cases, it is possible to obtain 
stationary states with finite energy as for the case of inelastic interactions. In 
such general framework, the corresponding initial value problem model is 



p{v,t) = at)Q{fJ){v,t) + g{f{t,v)), ^2.9) 
f{v,0) = fo{v), 

where the collision operator Q{f,f)iv,t) is as in (2.1) and Q{f(t,v)) models 
a 'heating source' due to different phenomena. The term ({t) may represent a 
mean field approximation that allows from proper time rescaling. See and 



15l | for several examples for these type of models and additional references. 



Following the work initiated in [15| and [IJ] on Non-Equilibrium Stationary 
States (NESS), our computational approach we shall present several compu- 
tational simulations of non-conservative models for either elastic or inelastic 
collisions associated to (2.9) of the Boltzmann equation with 'heating' sources. 
In all the cases we have addressed one can see that stationary states with finite 
energy are admissible, but they may not be Maxwellian distributions. Of this 
type of model we show computational output for three cases. First one is the 



pure diffusion thermal bath due to a randomly heated background [58|; |49|; |34 
in which case 

g,{f)=^MAf, (2.10) 
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where /i > is a constant. The second example relates to self-similar solutions 
of equation (2.9) for Q{f) = |45|;l25|, but dynamically rescaled by 

fM-^^n^MM). f = ^. (2.11) 



where 



Vo(t) = (a + r]t)-\ tm = -ln(l + -t), a, r] > 0. (2.12) 

1] a 

Then, the equation for f{v,t) coincides (after omitting the tildes) with equa- 
tion (2.9),for 

g2{f) = -vdiY{vf), v>0- (2.13) 
Of particular interest of dynamical time-thermal speed rescaling is the case 
of colhsional kernels corresponding to Maxwell type of interactions. Since the 
second moment of the collisional integral is a linear function of the energy, the 
energy evolves exponentially with a rate proportional to the energy production 
rate, that is 

^£{t) = Xo£{t), or equivalents ^(t) = £(0)6^"*, (2.14) 

with Ao the energy production rate. Therefore the corresponding rescaled vari- 
ables and equations (2.11) and (2. 9), (2. 13) for the study of long time behavior 
of rescaled solutions are 

f{v,t) = S-Ht)f{-^) = (^(0)e^«*)~i/>(^(0)e^»*)-^), (2.15) 

and / satisfies the self-similar equation (2.9) 
^2'(/) = —^o^fx, where x = vS~'^{t) is the self-similar variable. (2.16) 

We note that it has been shown that these dynamically self-similar states are 
stable under very specific scaling for a large class of initial states [12 . 

The last source type we consider is given by a model, related to weakly coupled 
mixture modeling slowdown (cooling) process [3] given by an elastic model in 
the presence of a thermostat given by Maxwell type interactions of particles 
of mass m having the Maxwellian distribution 

^-^^^=(2;^^ ' 

with a constant reference background or thermostat temperature T (i.e. the 
average of / Mq-dv = 1 and / \v\'^M'j-dv = T). Define 



QL{f){v, t) = / BL{\ul^l)f(v, t)Mr{'w) - f{v, t)Mr{w)] dadw . 

(2.17) 
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Then the corresponding evolution equation for f{v,t) is given by 



= QifJ)iv,t) + QQLif)iv,t) 
fiv,0) = fo{v). (2.18) 

where Q{f,f), defined as in (2.1), is the classical collision integral for elastic 
interactions (i.e. /3 = 1), so it conserves density, momentum and energy. The 
second integral term in (2.18) is a linear collision integral which conserves just 
the density (but not momentum or energy) since 



u = V — w the relative velocity 

V = V -\ {\u\a — u), w =w [\u\a — u) . 

m + 1 m + 1 

The coupling constant 9 depends on the initial density, the coupling constants 
and on m. The collision kernel Bi of the linear part may not be the same as 
the one for the non-linear part of the collision integral, however we assume that 
the Grad cut-off assumption (2.3) is satisfied and that, in order to secure mass 
preservation, the corresponding differential cross section functions and hi, 
the non-linear and linear collision kernels respectively, satisfy the renormalized 
condition 

bN{^) + QbL{^)da = 1 + 6. (2.20) 

Sd-l \U\ \U\ 

This last model describes the evolution of binary interactions of two sets of 
particles, heavy and light, in a weakly coupled limit, where the heavy particles 
have reached equilibrium. The heavy particles set constitutes the background 
or thermostat for the second set of particles. It is the light particle distribu- 
tion that is modeled by (2.18). Indeed, Q{f, f) corresponds to all the collisions 
which the light particles have with each other and the second linear integral 
term corresponds to collisions between light and heavy particles at equilibrium 
given by a classical distribution Mr{v). In this binary 3- dimensional, mixture 
scenario, collisions are assumed to be isotropic, elastic and the interactions 
kernels of Maxwell type. 

In the particular case of equal mass (i.e. m = 1), the model is of particu- 
lar interest for the development of numerical schemes and simulations bench- 
marks. Even though the local interactions are reversible (elastic), it does not 
conserve the total energy. In such there exists a special set of explicit, 

in spectral space, self-similar solutions which are attractors for a large class 
of initial states. When considering the case of Maxwell type of interactions 
in three dimensions i.e. B{\u\,fi) = b{fj,) with a cooling background process 
corresponding to a time temperature transformation, T = T(t) such that 



T(t) — > as t — i> 0, the models have self similar asymptotics [IJ; [l2] for a 
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large class of initial states. Such long time asymptotics corresponding to dy- 
namically scaled solutions of (2.18), in the form of (2.16), yields interesting 
behavior in f{v,t) for large time, converging to states with power like decay 
tails in v. In particular, such solution f{v,t) of (2.18) will lose moments as 



time grows, even if the initial state has all moments bounded, (see jl4j ; Il2l | for 
the analytical proofs). 



2.2 Collision Integral Representation 



One of the pivotal points in the derivation of the spectral numerical method 
for the computation of the non-linear Boltzmann equation lays in the repre- 
sentation of the collision integral in Fourier space by means of the weak form. 
Since for a suitably regular test function ip{v), the weak form of the collision 
integral takes the form (suppressing the time dependence in f) 



/ Q{f,mv)dv = 




f (v) f {w)B{\u\, fi)[ip{v')—'ilj{v)]dadwdv , 



(2.21) 

with v' = v + ^{\u\a — u). In particular, taking iljiv) = e (-\/27r)'^, where ( 
is the Fourier variable, we get the Fourier Transform of the collision integral 
through its weak form as 




']dadwdv . (2.22) 



We will use?= JF(.) - the Fourier transform and for the classical inverse 
Fourier transform. Plugging in the definitions of collision kernel B{\u\,fi) = 
C\{a)\u\'^ (which in the case of isotropic collisions would just be the Variable 
Hard Potential collision kernel) and v' 



= 77^ / . ^) / . ^^^)^^^ - u)e-^^-^dvdu 



GxAu,0[f{v)f{v-u)] du, (2.23) 



where 
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\u 



g5d-l 
A 



•Acres'*-! 



1^2 



(2.24) 



Note that (2.24) is valid for both isotropic and anisotropic interactions. For the 
former type, a simphfication ensues due to the fact the C\{a) is independent 
of (T G S'^-^: 



Gx,p{uX) = CxUd-2 \u\ 



,4 C.«sinc(^) - 1 



(2.25) 



Thus, it is seen that the dependence on a i.e. the integration over the unit 
sphere S"^~^ is completely done independently and there is actually a closed 
form expression for this integration, given by (2.25) in the case of isotropic 
collisions. In the case of anisotropic collisions, the dependence of C\ on a is 
again isolated into a separate integral over the unit sphere S'^~^ as given in 
(2.24). The above expression can be transformed for elastic collisions (3 = 1 
into a form suggested by Rjasanow and Ibragimov in their paper ^]. The 
corresponding expression for anisotropic collisions is given by (2.24). 

Further simplification of (2.23) is possible by observing that the Fourier trans- 
form inside the integral can be written in terms of the Fourier transform of 
f{v) since it can also be written as a convolution of the Fourier transforms. 

Let h{v) = fiy — u) 



1 



GxA^X) 



f{C-0f{0GxA^,0d^ 



fiC-OfiOe-'^-^'d^du 



(2.26) 



where GxA^X) = Lm'^ GxA'^X)e~'^"'du. Let u = re, e e S'^-\r e R For 
d = 3, ' 



GxA^X)= I I r^G{reX)e''^^-^dedT 

= IQn^Cx I r^+2[sinc(^^)sinc(r|e - -C|) - sinc(r|e|)]rfr . 

J r 2 2 



13 



Since the domain of computation is restricted to Q„ = [— L, L)^, u e [— 2L, 2L)^ =^ 
r e [0, 2 



JO 



■2^/3L 




A point worth noting is that the above formulation (2.26) results in 0{N'^'^) 
number of operations, where N is the number of discretizations in each velocity 
direction. Also, exploiting the symmetric nature in particular cases of the 
collision kernel one can reduce the number of operations to 0{NHogN). 



3 Numerical Method 

3. 1 Discretization of the Collision Integral 

Coming to the discretization of the velocity space, it is assumed that the two 
interacting velocities and the corresponding relative velocity 



where the velocity domain L is chosen such that u = v — w & LY 
through an assumption that supp{f) G [— L, L^. For a sufficiently large L, 
the computed distribution will not lose mass, since the initial momentum is 
conserved (there is no convection in space homogeneous problems), and is 
renormalized to zero mean velocity. We assume a uniform grid in the velocity 
space and in the fourier space with hy and /i^ as the respective grid element 
sizes, and /i^ are chosen such that hyh(^ — where N — number of 
discretizations of v and ( in each direction. 

3.2 Time Discretization 

To compute the actual particle distribution function, one needs to use an ap- 
proximation to the time derivative of /. For this, a second-order Runge-Kutta 
scheme or a Euler forward step method were used. Since a non-dimensional 
Boltzmann equation is computed, for numerical computations the value of 
time step dt is chosen such that it corresponds in its dimensional form to 
0.1 times the time between consecutive collisions (which depends on the col- 
lision frequency). During the standard process of non-dimensionalization of 
the Boltzmann Equation, such a reference quantity (time between collisions) 




(3.1) 
(3.2) 
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comes up. With time discretizations taken as = ndt, the discrete version of 
the Runge-Kutta scheme is given by 



(3.3) 



The corresponding Forward Euler scheme with smaller time step is given by 



3.3 Conservation Properties - Lagrange Multipliers 

Since the calculation of Qx,f3{f, f ){v) involves computing Fourier Transforms 
with respect to v, we extensively use Fast Fourier Transform. Note that the 
total number of operations in computing the collision integral reduces to the 
order of 3N'^Hog{N) +0{N^^) for (2.23) and 0{N^^) for (2.26). Observe that, 
choosing 1/2 < /3 < 1, the proposed scheme works for both clastic and inelas- 
tic collisions. As a note, the method proposed in the current work can also be 
extended to lower dimensions in velocity space. 

In the current work, due to the discretizations and the use of Fourier Trans- 
form, the accuracy of the proposed method relies heavily on the size of the 
grid and the number of points taken in each velocity/ Fourier space directions. 
Because of this it is seen that the computed Qx^f, f]{v) does not really con- 
serve quantities it is supposed to i.e. p, m, e for elastic collisions, p for Linear 
Boltzmann Integral and p, m for inelastic collisions. Even though the differ- 
ence between the computed (discretized) collision integral and the continuous 
one is not great, it is nevertheless essential that this issue be resolved. To rem- 
edy this, a simple constrained Lagrange multiplier method is employed where 
the constraints are the required conservation properties. Let M — A^'*, the to- 
tal number of discretizations of the velocity space. Assume that the classical 
Boltzmann collision operator is being computed. So p, m = (ml,m2,m3) and 
e are conserved. Let ujj be the integration weights where j — 1, 2, M. Let 



fiv^)^r{v^)+dtQ{r,n 



(3.4) 




be the distribution vector at the computed time step and 
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be the corrected distribution vector with the required moments conserved. Let 



C, 



(d+2)xM 



ViUj 



\ 



and 



''(d+2)xi - yP ml m2 m3 e 

be the vector of conserved quantities. Using the above vectors, the conservation 
can be written as a constrained optimization problem: 



* < 



11/ - /Hi ^ mm 

Cf = a; C e R'^+2x^, / e ^ g ^^+2 



To solve (*), one can employ the Lagrange multiplier method. Let A e 

be the Lagrange multiplier vector. Then the scalar objective function to be 

optimized is given by 



M 



Lif,X) = Y.\fj-fj\' + ^''{Cf-a). 



(3.5) 



Equation (3.5) can actually be solved exphcitly for the corrected distribution 
value and the resulting equation of correction be implemented numerically in 
the code. Taking the derivative of L(/, A) with respect to fj,j — 1, M and 
Aj, i = 1, ...,d + 2 i.e. gradients of L, 



^ = 0;j = l,...,M 
dfj 



And 



(3.6) 



dL 



Cf^a 

i.e. retrieves the constraints. 
Solving for A, 



CC^X = 2(a - Cf) . 



(3.7) 



(3.^ 
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Now CC"^ is symmetric (CC^)^ = CC^ and because C is the integration 
matrix, CC^ is positive definite. By linear algebra, the inverse of CC^ exists. 
In particular one can compute the value of A by 

\^2{CC^)-\a-Cf). (3.9) 

Substituting A into (3.6), 

f = f + C^^CC^)-\a-C]). (3.10) 

Using equation for forward Euler scheme (3.4), the complete scheme is given 

by (r(^^^■) = /;)VJ: 
/,=/^+rf^w^/^) 

/n+l^/. + C'^(CC^)-l(a_C/). 

So, 

= /r + dtQUJ, /7) + C^{CC^)-\a - Cf) 

= /; + dtQif^, /;) + C'iCC^r'ia - a - dtCQif^, /;)) 

= /; + dtQiff, /;) - dtc^{cc^)-'CQ{ff, /;) 
= /; + dt[i - c^{cc''r'c]Q{f^, /;) , 

(3.12) 

with I - X identity matrix. Letting Ajv(C) = I - C^{CC^)-^C with I - 
N X N identity matrix, one obtains 

= /; + dtAA,(c)g(/;, /;) , (3.13) 

where we expect the required observables are conserved and the solution ap- 
proaches a stationary state, since lim„_,.oo ||A7v(C) (5(/", /")||oo — . 

Identity (3.13) summarizes the whole conservation process. As described pre- 
viously, setting the conservation properties as constraints to a Lagrange mul- 
tiplier optimization problem ensures that the required observables are con- 
served. Also, the optimization method can be extended to have the distribu- 
tion function satisfy more (higher order) moments from (2.8). In this case, 
a{t) will include entries of ninit) from (5.1). 

We point out that for the linear Boltzmann collision operator used in the 
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mixture problem conserves density and not momentum(unless one computes 
isotropic solutions) and energy. For this problem, the constraint would just 
be the density equation. For inelastic collisions, density and momentum are 
conserved and for this case the constraint would be the energy and momen- 
tum equations. And for the elastic Boltzmann operator, all three quantities 
(density, momentum and energy) are conserved and thus they become the 
constraints for the optimization problem. The behavior of the conservation 
correction for Pseudo-Maxwell Potentials for Elastic collisions will be numer- 
ically studied in the numerical results section. This approach of using La- 
grangian constraints in order to secure moment preservation differs from the 



one proposed in [27|], [28|] for spectral solvers. 



4 Self-Similar asymptotics for a general elastic or inelastic BTE of 
Maxwell type or the cold thermostat problem - power law tails 



As mentioned in introduction, a new interesting benchmark problem for our 
scheme is that of a dynamically scaled solutions or self-similar asymptotics. 
More precisely, we present simulation where the computed solution in properly 
scaled time approaches a self similar solution. This is of interest because of the 
power tail behavior i.e. higher order moments of the computed solution are 
bounded. For the completeness of this presentation, the analytical description 
of such asymptotics is given in the following two sub sections. 



4-1 Self-Similar Solution for a non-negative Thermostat Temperature 



We consider the Maxwell type equation from (2.18) related to a space homo- 
geneous model for a weakly coupled mixture modeling slowdown process. The 
content of this section is dealt in detail in [l3] for a particular choice of zero 
background temperature (cold thermostat). For the sake of brevity, we refer 
to [14| for details. However, a slightly more general form of the self-similar 
solution for non zero background temperature is derived here from the zero 
background temperature solution. Without loss of generality for our numerical 
test, we assume the differential cross sections bi for collision kernel of the lin- 
ear and b^, the corresponding one for the nonlinear part, are the same, both 
denoted by b{^), satisfying the Grad cut-off conditions (2.3). In particular, 
condition (2.20) is automatically satisfied. 

In [l3], Fourier transform of the isotropic self-similar solution associated to 
problem in (2.18) will take the form: 

(f){x,t) = ^{xe-^"*) = l-a{xe-''y, as xe"''* ^ 0, with p<l, (4.1) 
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where x = |CP/2 and /i and 9 are related by 



2 (3p+l)(2-p) 
^=3^ ® = • 

Note that p = 1 corresponds to initial states with finite energy. It was shown in 



ij] for T = (i.e. cold thermostat), the Fourier transform of the self-similar, 



isotropic solutions of (2.18) is given by 



and its corresponding inverse Fourier transform, both for p = 1, /i = | and 
6 = I (as computed in [l4!]) is given by 

fn\v\,t)=e'Fo{\v\e'/') with Fo{\v\) = - / ds. (4.3) 

TT Jo (1 + s^y (27rs2)2 

Remark: It is interesting to observe that, as computed originally in for 
p = ^ or p = ^ in (4.2) yields 9 = 0, and one can construct explicit solutions 
to the elastic BTE with infinite initial energy. It is clear now that in order 
to have self-similar explicit solutions with finite energy one needs to have this 
weakly couple mixture model for slowdown processes, or bluntly speaking the 
linear coUisional term added to the elastic energy conservative operator. 

Finally, in order to recover the self-similar solution for the original equilib- 
rium positive temperature T (i.e. hot thermostat case) for the linear coUisional 
term, we denote, including time dependence for convenience. 



(f>o{x,t) = 4>{x,t)Thermostat=0 and (l)r{.X,t) = (f){x , t)Thermostat=T 

SO that (f)r{x,t) = 0o(a;, t)e~^^' . (4.4) 

Note that the solution constructed in (4.2) is actually (f)Q{x,t). Then the self- 
similar solution for non zero background temperature, denoted by (f)r{x,t) 
satisfies 



Mk, t) = - r e-\'\'^-'"'^^'''—^e-\'\'^/'ds 

Tl Jo (1 + S^)^ 

= - r e-^'^'^^~''''^^'^'^/'^-^ds. (4.5) 

71 Jo (1 + S^)"' 

In particular, let T = e~^*/'^as^ + T then, taking the inverse Fourier Trans- 
form, we obtain the corresponding self-similar state, according to (2.15), in 
probability space 
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f^^i\v\,t) = e'Fri\v\e'/') with Fri\v\) = - / ds. (4.6) 

vr Jo (1 + s-'y (27rT)2 

Then, letting t ^ oo, since T = T + as^e^ T, yields 



71^ (27rT)2 Jo (1 + 5"=)^ 

4 1 2 f s \ 
since — / 7:; ^c^s = — + arctan(s) I?" = 1. (4.8) 

So, the self-similar particle distribution /-f*(f , t) approaches a rescaled Maxwellian 
distribution with the background temperature T, that is according to (2.15), 



/^^(|^|,t) = e*Fr(|^|e*/^)^-4-^e-(H^^^*^^)/^^+*, as t ^ oo .(4.9) 

(27rT) : 



Remark: As pointed out in the previous remark, such asymptotic behavior, 
for finite initial energy, is due to the balance of the binary term and the linear 
coUisional term in (2.18). 



In addition, very interesting behavior is seen on Fr(|f |) as T — (cold 
thermostat problem), where the particle distribution approaches a distribu- 
tion with power-like tails (i.e. a power law decay for large values of |f |) and 
an integral singularity at the origin. Indeed, in 1^ an asymptotic behavior is 
derived for -fo(l^l) from (4.3), for large and small values of \v\, leading to 



F{\v\) = 2{-)'/'-^^[l + 0{^)], for |t;|^oo, 
vr l^;!" |f| 

F{\v\) = ^7^[l + 2\v\Hn{\v\) + 0{\v\% for |'^;| ^ 0. (4.10) 

In particular the self-similar particle distribution function F(|f|), v G M^, 
behaves like as \v\ oo, and as ■j;;^^ as |f | — 0, which indicates a very 
anomalous, non-equilibrium behavior as function of velocity; but, nevertheless, 
remains with finite mass and kinetic temperature. This asymptotic effect can 
be described as an overpopulated (with respect to Maxwellian), large energy 
tails and infinitely many particles at zero energy. This interesting, unusual 
behavior is observed in problems of soft condensed matter js^] • 



We shall see, then in the following section, that our solver captures these 
states with spectral accuracy and consequently the self similar solutions are 
attractors for a large class of initial states. These numerical tests are a cru- 
cial aspect of the spectral Lagrangian deterministic solver used to simulate 
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this type of non-equilibrium phenomena, where all these explicit formulas for 
our probability distributions allow us to carefully benchmark the proposed 
numerical scheme. 



4-2 Self-Similar asymptotics for a general problem 



The self-similar nature of the solutions |) for a general class of problems, 
for a wide range of values for the parameters /?, p, fi and 6, was addressed 
in [lil with much detail. Three different behaviors have been clearly explained. 
Of particular interest for our present numerical study are the mixture problem 
with a cold background and the inelastic Boltzmann cases. Interested readers 



are referred to 12 



For the purpose of our presentation, let = J-'[f] be the Fourier transform of 
the probability distribution function satisfying the initial value problem (2.1)- 
(2.2) or (2.9). Let's denote by r(0) = J^[Q~^{f, f)] the Fourier transform of 
the gain part of the collisional term associated with the initial value problem. 
It was shown in that the operator T{(f)), defined over the Banach space of con- 
tinuous bounded functions with the L°°-norm (i.e. the space of characteristic 
functions, that is the space of Fourier transforms of probability distributions), 
satisfies the following three properties [12i] : 

1 - r(0) preserves the unit ball in the Banach space. 

2 - r(0) is L-Lipschitz operator, i.e. there exists a bounded linear operator L 

in the Banach space, such that 

|r(«i)-r(M2)|(x,t) < L(|Mi-M2|(x,t)), V \\Ui\\<l;t = l,2. 

3 - r(0) is invariant under transformations (dilations) 

e^^T(u) = T(e^^u) , V = x-^ , e^^u(x) = u(xe^), r G M+ . (4.11) 

ox 

In the particular case of the initial value problem associated to Boltzmann 
type of equations for Maxwell type of interactions, the bounded linear opera- 
tor that satisfies property 2, is the one that linearizes the Fourier transform 
of the gain operator about the state u = 1. 

Next, let x^ be the eigenfunction corresponding to the eigenvalue X{p) of the 
linear operator L associated to F, i.e. L{x^) = X{p)x^. 

Define the spectral function associated to F given by /i(p) = ^^^^ defined 
for p > 0. It was shown in [l2| that A'.(0+) = +oo (i.e. p = is a vertical 



asymptote) and that for the problems associated to the initial value problems 
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(2.1)-(2.2) or (2.9), there exists a unique minimum for fi{p) localized at po > 1, 
and that fi{p) — 0~ as p ^ +00. 

Then, the existence of self-similar states and convergence of the solution to the 
initial value problem to such self-similar distribution function was described 



in [12] in the following four statements: 



(i) [13 / - Lemma 5.1 (existence): There exists a unique isotropic solution f{\v\,t) 
to the initial value problem (2.1)-(2.2) or (2.9) for Maxwell type interac- 
tions, in the class of probability measures, satisfying /(|f|,0) = /o(|f|) > 
0, J^d fo{\v\)dv = 1 such that for the Fourier transform problem x = Uq = 
nfo{\v\)] = l + 0{x),asx^O, 



(ii) Self similar states - Theorem 7.2: f{\v\,t) has self-similar asymptotics in the 
following sense: 

Taking the Fourier transform of the initial state to satisfy 

Uo + fi{p) xP Uq = T{uo) + 0{xP+^), such that p + e < po , (4.12) 

(i.e. fi{p);fi'{p) < 0). Then, there exists a unique, non-negative, self-similar 
solution 

ri\v\,t) = e-t^(^')*Fp(|t;|e-^^(^')*) , 
with J^{Fp{\v\)) = w{x),x = 72 s.t. fi{p)xPw'{x) + w{x) = T{w). 



(iii) Self similar asymptotics - Section 9 and Theorem 11.1 in \12] : There exists 
a unique (in the class of probability measures) solution f{\v\,t) satisfying 
/(|f 1,0) = /o(|f I) > 0, with /^^ f{\v\)dv = 1, such that, for x = and 

T[h{\v\)] = l-ax^ + 0{x'P^'),x^Q,Q<p<l withp + e<po • 



Then, for any given < p < 1, there exists a unique non- negative self-similar 
solution /i?)(|w|,t) = e-i^(p)*Fp(|w|e-^A'(p)*) such that 

f{\vlt) ^,^00 e-i'^(^)*Fp(|t;|e-5^(^)*) . (4.13) 

or equivalently 

elMp)*/(|^|eiMp)t^i) ^^(1^1) ^ (4.14) 

where ii{p) is the value of spectral function associated to the linear bounded 
operator L as described above. 
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Fig. 1. Spectral Function for a general homogeneous Boltzmann collisional 
problem of Maxwell type 

(iv) Power tail behavior of the asymptotic hmit: If fi{p) < 0, then the self- 
similar limiting function Fp(|f |) does not have finite moments of all orders. 
In addition, if < p < 1 then all moments of order less than p are bounded; 
i.e. rUg = J^d Fp{\v\)\v\'^'^dv < C!c;0 < q < p. However, if p = 1 (finite energy- 
case) then, the boundedness of moments of any order larger than 1, depend 
on the conjugate value of by the spectral function That means 

rriq < oo only for < g < , where p^ > po > 1 is the unique maximal root 
of the equation = 

Remark 1: When p = 1, is the energy dissipation rate, and £{t) = e^'-^^* 
the kinetic energy evolution function. So, £{tY^^f{v£{t),t) Fi{\v\). 

Remark 2: We point out that condition (4.12) on the initial state is eas- 
ily satisfied by taking a sufficiently concentrated Maxwellian distribution as 
shown in fl^ l. and as done for our simulations in the next section. 

However, rescaling with a different rate, it is not possible to pick up the non- 
trivial limiting state Z^**, since 

/(|^|e>,t)^,^^e->5o(|^|); r/>/x(l), (4.15) 

and 

fi\v\e-^''\ t) ^t-.o. 0; /i(p™„) <fiil + 5)<r]< /i(l) . (4.16) 
These results are also true for any p < I. For the general space homogeneous 
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(elastic or inelastic) Boltzmann model of Maxwell type or the corresponding 
mixture problem, the spectral function is given in Figure 1. 



5 Numerical Results 

We benchmark the new proposed numerical method to compute several ex- 
amples of 3 — D in velocity and time for initial value problems associated with 
non-conservative models where some analysis is available, as are exact mo- 
ment formulas for Maxwell type of interactions as well as qualitative analysis 
for solutions of VHS models. We shall plot our numerical results versus the 
exact available solutions in several cases. Because all the computed problems 
converge to an isotropic long time state, we choose to plot the distribution 
function in only one direction, which is chosen to be the one with the ini- 
tial anisotropics. All examples considered in this manuscript are assumed to 
have isotropic, VHS collision kernels, i.e. differential cross section indepen- 
dent from scattering angle. We simulate the homogeneous problem associated 
to the following problems for different choices of the parameters j3 and A, and 
the Jacobian and heating force term Q{f ). 

5. 1 Maxwell type of Elastic Collisions 

Consider the initial value problem (2.1), (2.2), with B{\u\,ij) — In 
(2.1), (2.2), the value of the parameters are P = 1, Jp = 1 and A = with 
the pre-collision velocities defined from (2.2). In this case, for a general initial 
state with finite mass, mean and kinetic energy, there is no exact expression 
for the evolving distribution function. However there are exact expressions 
for all the statistical moments (observables) . Thus, the numerical method is 
compared with the known analytical moments for different discretizations in 
the velocity space. 

The initial states we take are convex combinations of two shifted Maxwellian 
distributions. So consider the following case of initial states with unit mass 
Jjjs fQ{v)dv — 1 given by convex combinations of shifted Maxwellians 

f{v, 0) = /o(t) = ^MtAv - Vi) + (1 - 7)Mt,{v - V2); with ^ 7 ^ 1 

-\y-V\'^ 

where Mt{v — V)— (27rr)a/2 ^ '^'^^ • Then, taking 7 = 0.5 and mean fields for 
the initial state determined by 

yi = [-2,2,0r, 1/2 = [2,0, of; = 1 ,T2 = 1 , 
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enables the first five moment equations corresponding to the collision invari- 
ants to be computed from those of the initial state. All higher order moments 
are computed using the classical moments recursion formulas for Maxwell type 
of interactions (2.8). In particular, it is possible to obtain the exact evolution 
of moments as functions of time. Thus 



p{t) = po = 1 and V{t) = = [0, 1, 0]^ . 

By a corresponding moment calculation as in (2.8), the complete evolution of 
the second moment tensor (2.7) is given by 



M{t) 



and the energy flow (2.7) 



r{t) = 2 
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-t/3 



^l2^ 



(e-*/2 - e 



and the kinetic temperature is conserved, so 

m = To = \. 



(5.1) 



The above moments along with their numerical approximations for different 
discretizations in velocity space are plotted in Figures 5.1. In Figure 3, the 
evolution of the computed distribution function into a Maxwellian is plotted 
for = 40. In order to check the conservation accuracy of the method, let 
- unconserved distribution given as input to the conservation routine and fc - 
conserved distribution resulting from the conservation routine. With a convex 
combination of two Gaussians as input, the numerical method is allowed to 
run and ||/c — /m||oo is plotted for all times for different values of in figure 
4. As expected, for t approaching the final time, the largest value of A^ gives 
the smallest conservation correction. 



5.2 Maxwell type of Elastic collisions - Bobylev-Krook-Wu (BKW) Solution 



An explicit solution to the initial value problem (2.1) for elastic. Maxwell 
type of interactions (/? = 1, A = 0) was derived in and independently in 
431 1 for initial states that have at least 2 + (5-moments bounded. It is not of 
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Fig. 2. Maxwell type of Elastic collisions: Momentum Flow Afn , M12 , M22 , M33 , 
Energy Flow ri , r2 

self-similar type, but it can be shown to converge to a Maxwellian distribution. 
This solution takes the form 
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Fig. 3. Maxwell type of Elastic collisions: Evolution of the Distribution function 





Fig. 4. Maxwell type Elastic Collisions: Conservation Correction for Elastic Colli- 
sions 

where K = 1 — e^*/^ and rj =initial distribution temperature. It is interesting 
that it is negative for small values of t. So in order to obtained a physically 
meaning probability distribution, / must be non-negative. This is indeed the 
case for any K ^ ^ or t ^ to = 6/n(|) ~ 5.498. In order to test the accuracy 
of our solver, set the initial distribution function to be the BKW solution, 
the numerical approximation to the BKW solution and the exact solution are 
plotted for different values of at various time steps in Figure 5. 




N = 8 
N = 16 
N = 32 
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V_x V_x 

(A^ = 24) (A^ = 32) 

Fig. 5. BKW, p,E{t) conserved 

5.3 Hard-Sphere Elastic Collisions 

In (2.1), (2.2), we have /5 = 1, Jg = 1 and A = 1 with the post-colhsion veloci- 
ties defined from (2.1). Unhke Maxwell type of interactions, there is no explicit 
expression for the moment equations and neither is there any explicit solution 
expression as in the BKW solution scenario. For Hard Sphere isotropic colli- 
sions, the expected behavior of the moments is similar to that of the Maxwell 
type of interactions case except that in this case, the moments somewhat 
evolve to the equilibrium a bit faster than in the former case i.e. figure 6. Also 
plotted is the time evolution of the distribution function starting from the 
convex combination of Maxwellians as described in a previous subsection in 
Figure 7. 



5.4 Inelastic Collisions 



This is the scenario wherein the utility of the proposed method is clearly seen. 
No other deterministic method can compute the distribution function in the 
case of inelastic collisions (isotropic), but the current method computed this 
3~D evolution without much complication and with the exactly same number 
of operations as used in an elastic collision case. This model works for all sorts 
of Variable Hard Potential interactions. Consider the special case of Maxwell 
(A = 0) type of inelastic (/? 7^ 1) collisions in a space homogeneous Boltzmann 
Equation in (2.1), (2.2). Let (f){v) = |f p be a smooth enough test function. 
Using the weak form of the Boltzmann equation with such a test function one 
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Fig. 6. Hard Sphere, Elastic: Momentum Flow Mn, M12, M22, M33, Energy Flow 
ri,r2 



can obtain the ODE governing the evolution of the kinetic energy K{t) 



(5.3) 
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Fig. 7. Hard-Sphere, Elastic: Evolution of the Distribution function, N = 32 

where V - conserved (constant) bulk velocity of the distribution function. This 
gives the following solution for the kinetic energy as computed in (2.8) 

Kit) = K(0)e-^(i-^)* + ^(1 - e-'^^^-'^)*) , (5.4) 

where K{0) = kinetic energy at time t = 0. As we have an explicit expression 
for the kinetic energy evolving in time, this analytical moment can be com- 
pared with its numerical approximation for accuracy and the corresponding 
graph is given in Figure 8. Also the general evolution of the distribution in 
an inelastic collision environment is also shown in Figure 8. In the conser- 
vation routine (constrained Lagrange multiplier method), energy is not used 
as a constraint and just density and momentum equations are used for con- 
straints. Figure 8 shows the numerical accuracy of the method even though 
the energy (plotted quantity) is not being conserved as part of the constrained 
optimization method. 



5.5 Inelastic Collisions with Diffusion Term 

Here we simulate, the equations (2.9), (2.10). Here we simulate a model cor- 
responding to inelastic interactions in a randomly excited heat bath with 
constant temperature rj. The evolution equation for kinetic temperature as 
a function of time is given by: 

dT 1 — /■ r r 

= 2r^-C—— (l-i^)B{\ulfj.)\u\'f{v)fiw)dadwdv, 
at 2A Jvm^ Jwm^ Jaes^ 

(5.5) 
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Fig. 8. Inelastic: Kinetic Energy (left) & f{v, t) (right) 

which, in the case of inelastic Maxwell type of interactions according to (2.8), 
(5.5) becomes 

^ = 2r^ - (nC,{l - e')T . (5.6) 

The above equation gives a closed form expression for the time evolution of 
the kinetic temperature and can be expressed as follows: 



where 



To = - / \v\'f{v)dv and - 



3 A.GR3 ' ' ^ ^ ^ ~ CttCo{1 - e2 



As it can be seen from the expression for T, in the absence of the diffusion 
term (i.e. rj = 0) and for e ^ 1 (inelastic collisions), the kinetic temperature of 
the distribution function decays like an exponential just like in the previous 
section. So, the presence of the diffusion term pushes the temperature to an 
equilibrium value of > even in the case of inelastic collisions. Also 

note that if the interactions were elastic and the diffusion coefficient positive 
then, = +00, so there would be no equilibrium states with finite kinetic 



temperature. These properties were shown in [3J| and similar time asymptotic 
behavior is expected in the case of hard-sphere interactions where > 
is shown to exist. However, the time evolution of the kinetic temperature is a 
non-local integral (5.5) does not satisfy a close ODE form (5.6). The proposed 
numerical method for the calculation of the collision integral is tested for 
these two cases. We compared with the analytical expression (5.7) for different 
initial data, the corresponding computed kinetic temperatures for Maxwell 
type interactions in Figure 9. The asymptotic behavior is observed in the case 
of hard-sphere interactions in Figure 10. The conservation properties for this 
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Fig. 9. Maxwell type of Inelastic collisions, Diffusion Term for N = 16 




Fig. 10. Hard-Sphere, Inelastic Collisions, Diffusion Term, < Tq for N = 16 

case of inelastic collisions with a diffusion term are set exactly like in the 
previous subsection (inelastic collisions without the diffusion term). 



5. 6 Maxwell type of Elastic Collisions - Slow down process problem 



Next, consider (2.18) with /5 = 1, Jg = 1 and B{\u\,n) = ^ i.e. isotropic 
collisions. The second term is the linear collision integral which conserves 
only density and the the first term is the classical collision integral from 
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(2.18) conserving density, momentum and energy. M{v) in (2.18) refers to 

the Maxwellian defined by Mr{v) = e (^t^) ^^-^^372, with T the constant ther- 
mostat temperature. In particular, any initial distribution function converges 
to the background distribution Mr- Such behavior is well captured by the 
numerical method. Indeed, Fig. 11 corresponds to an initial state of a convex 
combination of two Maxwellians. In addition, from (4.6): 



\/(2) /-oo 1 p-blV2f 



which is the finite energy solution for p = 1, a = 1, = |, 6* = | in (4.2), i.e. 
p = 1 in (4.13) and (4.14). As t — >■ 00, the time rescaled numerical distribution 
is compared with the analytical solution /|-* for a positive background tem- 
perature T and it converges to a Maxwellian Mq-- From Figure 11, it can be 
seen that the numerical method is quite accurate and the computed distribu- 
tion is in very good agreement with the analytical self-similar distribution fj-^ 
from (4.6). Similar agreement is observed for different constant values of T 
approaching (Figure 11). The interesting asymptotics (4.10) corresponding 
to power-like tails and infinitely many particles at zero energies occur only 
when T = as shown in (4.10) and (4.10). 

Since letting T = in the scheme created an instability, we proposed the 
following new methodology to counter this effect. We let, instead, T = (e'""^ 
ensuring that the thermostat temperature vanishes for large time and set 



T = Ce""* + as'^e^ , (5.8) 

where the role of a is very important and a proper choice needs to be made. 
In our simulations, we take ( = 0.25 and the values of a need to be chosen 
exactly as a = /i(l) = 2/3, the energy dissipation rate as described in section 
4.2 to recover the asymptotics as in (4.10). 

Remark: Due to the exponential time rescaling of Fourier modes, our proce- 
dure to compute self-similar solutions in free space may also be viewed as a 
non-uniform grid of Fourier modes that are distributed according to the con- 
tinuum spectrum of the associated problem. This choice plays the equivalent 
role to the corresponding spectral approximation of the free space problem 
of the heat kernel, that is, the Green's function for the heat equation, which 
happens to be a similarity solution as well, due to the linearity of the problem 
in this case. In particular, we expect optimal algorithm complexity using such 
non-equispaced Fast Fourier Transform, as obtained by Greengard and Lin 



3a] for spectral approximation of the free space heat kernel. This problem 



will be addressed in a forthcoming paper. 
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Fig. 11. Maxwell type collisions, Slow down process with B = 4/3, /i = 2/3, N = 24 



The following plots elucidate the fact that power-like tails are achieved asymp- 
totically with a decaying T. For a decaying background temperature as in 
(5.8), Figure 12 shows evolution of a convex combination of Maxwellians to a 
self-similar (blow up for zero energies and power-like for high energies) behav- 
ior. Figure 13 plots the computed distribution along with a Maxwellian with 
temperature of the computed solution. This illustrates that the computed self- 
similar solution is largely deviated from a Maxwellian equilibrium. In order 
to better capture the power-like effect using this numerical method, we set 
T = Ce~^*/^ = C^"'^*, see (5.8), where /i is related the spectral properties of 
the Fourier transformed equation as described in section 4.2 on the slow down 
process problem with fj, = fj,{l) the energy dissipation rate. Thus, as it was 
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Fig. 12. Slow down process: N = 32, T 
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Fig. 13. Computed distribution Vs. Maxwellian with temperature of the computed 
distribution 



computed in [IJ] and revised in section 4 of this paper, we know that for initial 
states with finite energy, p = 1 and the corresponding energy dissipation rate 
is fj,{l) = /i = 2/3 is positive. In particular p^, = 1.5 is the conjugate of p = 1 of 
the spectral curve nig in Theorem 4.1 part (i). In addition the rescaled proba- 
bility will converge to the moments of the self-similar state (4.13), (4.14), that 
is 

-'?*2/3 / f{v)\v\^''d 



and we know any moment m„ is unbounded for g > = 1.5. 
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We have plot in Fig. 14 the evolution of e ''^'^^'^ Jy^^3 f{v)\v\'^''dv for q = 
1, 1.3, 1.45, 1.5, 1.55, 1.7, 2.0, computed for different values of iV = 10, 14, 16, 
18, 22, 26. It can be seen that as time progresses (and as the thermostat tem- 
perature T decreases to 0) , the approximated numerically computed moments 
to niq, q > 1.5 start to blow up as predicted. The value q — 1.5 is the thresh- 
old value, as any moment mq>i,5{t) — > oo. The expected spectral accuracy, 
as the value of increases, improves the growth zone of such moments for 
larger final times. The reason for such effect is that since the truncation of 
Fourier modes that results in the truncation in velocity domain, makes the 
distribution function to take small negative values for large velocities con- 
tributing to numerical errors that may cause rUq to peak and then relax back. 
In particular, larger order moments of the computed self-similar asymptotics 
with the negative oscillating parts on large energy tails, result in the large 
negative moment values for the above mentioned values of N crating large 
negative errors. However it is noticed that the negative oscillation values of 
f{t, v) coincide with large velocity values used in getting its g-moments ap- 
proximating rriq, for q > 1.5, and that such error is reduced in time for larger 
number of Fourier modes. Finally we point out that a FFTW package has 
been used. We have noticed in our numerics that for the specific choice of 
values 7^ 6, 10, 14, 18, 22, 26, 6 + 4A;; A; = 0, 1, 2, 3, the approximating 
moments to mg(t) start to take negative values very quickly, as seen in Fig. 15 
for = 16 and 20, making the numerical solution inadmissible since analyt- 
ically rriqit) > 0,Vt. Such effect may be due to the particular choice of the 
FTTW solver. 



6 Conclusions and Future Work 

In conclusion, the presented numerical method works for elastic and inelastic 
variable hard potential interactions. This is first of its kind as no additional 
modification is required to compute for elastic and inelastic collisions. In com- 
parison with the known analytical results (moment equations for elastic BTE, 
BKW self-similar solution, attracting Bobylev-Cercignani-Gamba self-similar 
solutions for elastic collisions in a slow down process), the computed ones are 
found to be very close. The method employs a Fast Fourier Transform for faster 
evaluation of the collision integral. Even though the method is implemented 
for a uniform grid in velocity space, it can even be implemented for a non- 
uniform velocity grid. The only challenge in this case is computing the Fast 
Fourier Transform on such a non-uniform grid. There are available packages 
for this purpose, but such a non-uniform FFT can also be implemented using 
certain high degree polynomial interpolation and this possibility is currently 
being explored. The integration over the unit sphere is avoided completely and 
only a simple integration over a regular velocity grid is needed. Even though 
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= 10 



N = U 




a trapezoidal rule is used as an integration rule, other integration rules like a 
Gaussian quadrature can be used to get better accuracy. For time discretiza- 
tion, a simple second-order Runge Kutta scheme is used. The proposed method 
has a big advantage over other non-deterministic methods as the exact distri- 
bution function can actually be computed instead of just the averages. 

Implementations of this scheme for the space inhomogeneous case is currently 
developed by the authors by means of splitting algorithms in advection and 
collision components. Next step in this direction would be to implement the 
method for a practical 1 and 2 — D space inhomogeneous problems such shock 
tube phenomena for specular and diffusive boundary conditions, and Rayleigh- 
Benard instability or a Couette flow problem. 
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Fig. 15. mq{t) for T = e"^*/^ 

7 Acknowledgements 

The authors would like to thank Sergej Rjasanow for discussions about the 
conservation properties of the numerical method and for other comments. 
Both authors are partially supported under the NSF grant DMS-0507038. 
Support from the Institute of Computational Engineering and Sciences and 
the University of Texas Austin is also gratefully acknowledged. 



References 

[1] R. J. Alonso, I. M. Gamba, Propagation of and maxwellian weighted 
bounds for derivatives of solutions to the homogeneous elastic boltzmann 
equation, To appear in Journal of Mathematiques Pures et Appliquees. 

[2] G. A. Bird, Molecular Gas Dynamics, Clarendon Press, Oxford, 1994. 

[3] A. V. Bobylev, Exact solutions of the boltzmann equation, (Russian) 
Dokl. Akad. Nauk SSSR 225 (1975) 1296-1299. 

[4] A. V. Bobylev, Exact solutions of the nonlinear boltzmann equation and 
the theory of relaxation of a maxwellian gas. Translated from Teoretich- 
eskaya i Mathematicheskaya Fizika 60 (1984) 280 - 310. 

[5] A. V. Bobylev, The theory of the nonlinear spatially uniform boltzmann 
equation for maxwell molecules. Mathematical physics reviews; Soviet 
Sci. Rev. Sect. C Math. Phys. Rev. 7 (1988) 111-233. 

[6] A. V. Bobylev, J. A. Carrillo, I. M. Gamba, On some properties of ki- 
netic and hydrodynamic equations for inelastic interactions. Journal of 
Statistical Physics 98 (2000) 743-773. 



38 



[7] A. V. Bobylev, C. Cercignani, Discrete velocity models without nonphys- 

ical invariants, Journal of Statistical Physics 97 (1999) 677-686. 
[8] A. V. Bobylev, C. Cercignani, Exact eternal solutions of the boltzmann 

equation. Journal of Statistical Physics 106 (2002) 1019-1038. 
[9] A. V. Bobylev, C. Cercignani, The inverse laplace transform of some an- 
alytic functions with an application to the eternal solutions of the boltz- 
mann equation. Applied Mathematics Letters 15 (2002) 807-813(7). 
A. V. Bobylev, C. Cercignani, Moment equations for a granular material 
in a thermal bath. Journal of Statistical Physics 106 (2002) 547-567(21). 
A. V. Bobylev, C. Cercignani, Self-similar asymptotics for the boltzmann 
equation with inleastic and elastic interactions. Journal of Statistical 
Physics 110 (2003) 333-375. 

A. V. Bobylev, C. Cercignani, I. M. Gamba, On the self-similar asymp- 
totics for generalized non-linear kinetic maxwell models, to appear in 
Communication in Mathematical Physics. 
URL |http : //arxiv . org/abs/math-ph/ 0608035] 

A. V. Bobylev, C. Cercignani, G. Toscani, Proof of an asymptotic prop- 
erty of self-similar solutions of the boltzmann equation for granular ma- 
terials. Journal of Statistical Physics 111 (2003) 403-417. 
A. V. Bobylev, I. M. Gamba, Boltzmann equations for mixtures of 
maxwell gases: Exact solutions and power like tails. Journal of Statis- 
tical Physics 124 (2006) 497-516. 

A. V. Bobylev, I. M. Gamba, V. Panferov, Moment inequalities and high- 
energy tails for boltzmann equations with inelastic interactions, Journal 
of Statistical Physics 116 (2004) 1651-1682. 

A. V. Bobylev, M. Groppi, G. Spiga, Approximate solutions to the prob- 
lem of stationary shear flow of smooth granular materials, Eur. J. Mech. 
B Fluids 21 (2002) 91-103. 

A. V. Bobylev, S. Rjasanow, Difference scheme for the boltzmann equa- 
tion based on the fast fourier transform, European journal of mechanics. 

B. Fluids 16:22 (1997) 293-306. 

A. V. Bobylev, S. Rjasanow, Fast deterministic method of solving the 
boltzmann equation for hard spheres, European journal of mechanics. B, 
Fluids 18:55 (1999) 869-887. 

A. V. Bobylev, S. Rjasanow, Numerical solution of the boltzmann equa- 
tion using fully conservative difference scheme based on the fast fourier 
transform. Transport Theory Statist. Phys. 29 (2000) 289-310. 
J. E. Broadwell, Study of rarefied shear flow by the discrete velocity 
method, J. Fluid Mech. 19 (1964) 401-414. 

H. Cabannes, Global solution of the initial value problem for the discrete 
boltzmann equation. Arch. Mech. (Arch. Mech. Stos.) 30 (1978) 359-366. 

C. Cercignani, Recent developments in the mechanics of granular mate- 
rials, Fisica matematica e ingegneria delle strutture, Bologna: Pitagora 
Editrice (1995) 119-132. 

C. Cercignani, Shear flow of a granular material. Journal of Statistical 



39 



Physics 102 (2001) 1407-1415. 

[24] C. Cercignani, H. Cornille, Shock waves for a discrete velocity gas mix- 
ture, Journal of Statistical Physics 99 (2000) 115-140. 

[25] M. H. Ernst, R. Brito, Driven inelastic maxwell models with high energy 
tails, Phys. Rev. E 65 (4) (2002) 040301. 

[26] M. H. Ernst, R. Brito, Scaling solutions of inelastic boltzmann equations 
with over-populated high energy tails. Journal of Statistical Physics 109 

(2002) 407-432. 

[27] F. Filbet, C. Mouhot, L. Pareschi, Solving the boltzmann equation in 

nlogn, SIAM J. Sci. Comput. 28 (2006) 1029-1053. 
[28] F. Filbet, G. Russo, High order numerical methods for the space non 

homogeneous boltzmann equation. Journal of Computational Physics 186 

(2003) 457-480. 

URL citeseer . ist .psu. edu/f ilbetOShigh.html 

[29] F. Filbet, G. Russo, High order numerical methods for the space non 
homogeneous boltzmann equation. Journal of Computational Physics 186 
(2003) 457 - 480. 

[30] N. Fournier, S. Mischler, A boltzmann equation for elastic, inelastic and 

coalescing collisions. Journal de mathematiques pures et appliques 84 

(2005) 1173-1234. 
[31] M. Frigo, S. G. Johnson, Fast fourier transform of the west. 

URL www . f f tw . org 
[32] E. Gabetta, L. Pareschi, G. Toscani, Relaxation schemes for nonlinear 

kinetic equations, SIAM J. Numer. Anal. 34 (1997) 2168-2194. 
[33] I. M. Gamba, V. Panferov, C. Villani, Upper maxwellian bounds 

for the spatially homogeneous boltzmann equation. To appear in 

Arch.Rat.Mec.Anal. 

URL http : //arxiv.org/abs/math/07010811 

[34] I. M. Gamba, V. Panferov, C. Villani, On the boltzmann equation for 
diffusively excited granular media. Communications in Mathematical 
Physics 246 (2004) 503-541(39). 

[35] I. M. Gamba, S. Rjasanow, W. Wagner, Direct simulation of the uniformly 
heated granular boltzmann equation. Mathematical and Computer Mod- 
elling 42 (2005) 683-700. 

[36] I. M. Gamba, S. H. Tharkabhushanam, Convergence and error analysis 
of spectral-lagrange boltzmann solver. In preparation. 

[37] R. L. Greenblatt, J. L. Lebowitz, Product measure steady states of gen- 
eralized zero range processes, J. Phys. A 39 (2006) 1565-1573. 

[38] L. Greengard, P. Lin, Spectral approximation of the free-space heat ker- 
nel, Appl. Comput. Harmon. Anal. 9 (1) (2000) 83-97. 

[39] M. Herty, L. Pareschi, M. Seaid, Discrete-velocity models and relaxation 
schemes for traffic flows, SIAM J. Sci. Comput. 28 (2006) 1582-1596. 

[40] I. Ibragimov, S. Rjasanow, Numerical solution of the boltzmann equation 
on the uniform grid. Computing 69 (2002) 163-186. 

[41] R. lUner, On the derivation of the time-dependent equations of motion 



40 



for an ideal gas with discrete velocity distribution, J. de Mecanique 17 
(1978) 781-796. 

[42] S. Kawasliima, Global solution of the initial value problem for a discrete 
velocity model of the boltzmann equation, Proc. Japan Acad. Ser. A 
Math. Sci. 57 (1981) 19-24. 

[43] K. Max, W. T. Tsun, Formation of maxwellian tails. Physical Review 
Letters 36 (1976) 1107-1109. 

[44] L. Mieussens, Discrete-velocity models and numerical schemes for the 
boltzmann-bgk equation in plane and axisymmetric geometries. Journal 
of Computational Physics 162 (2000) 429-466. 

[45] J. M. Montanero, A. Santos, Computer simulation of uniformly heated 
granular fluids. Gran. Matt. 2 (2000) 53-64. 

[46] S. J. Moon, M. D. Shattuck, J. Swift, Velocity distributions and corre- 
lations in homogeneously heated granular media. Physical Review E 64 
(2001) 031303. 

[47] C. Mouhot, L. Pareschi, Fast algorithms for computing the boltzmann 

colhsion operator. Math. Comp. 75 (2006) 1833-1852. 
[48] K. Nanbu, Direct simulation scheme derived from the boltzmann equation 

i.monocomponent gases, J. Phys. Soc. Japan 52 (1983) 2042 - 2049. 
[49] T. V. Noije, M. Ernst, Velocity distributions in homogeneously cooling 

and heated granular fluids. Gran. Matt. 1:57(1998). 
[50] V. Panferov, S. Rjasanow, Deterministic approximation of the inelastic 

boltzmann equation. Unpublished manuscript. 
[51] L. Pareschi, B. Perthame, A fourier spectral method for homogenous 

boltzmann equations. Transport Theory Statist. Phys. 25 (2002) 369- 

382. 

[52] L. Pareschi, G. Russo, Numerical solution of the boltzmann equation. 

i. spectrally accurate approximation of the collision operator, SIAM J. 

Numerical Anal. (Online) 37 (2000) 1217-1245. 
[53] L. Pareschi, G. Toscani, Self-similarity and power-like tails in nonconser- 

vative kinetic models, J. Stat. Phys. 124 (2006) 747-779. 
[54] S. Rjasanow, W. Wagner, A stochastic weighted particle method for the 

boltzmann equation. Journal of Computational Physics (1996) 243-253. 
[55] S. Rjasanow, W. Wagner, Stochastic Numerics for the Boltzmann Equa- 
tion, Springer, Berlin, 2005. 
[56] C. Villani, Handbook of Fluid dynamics, chap. A Review of Mathematical 

Topics in Collisional Kinetic Theory, Elsevier, 2003, pp. 71-306. 

URL http : //www . umpa . ens-lyon . f r/~cvillani 
[57] W. Wagner, A convergence proof for bird's direct simulation monte carlo 

method for the boltzmann equation, Journal of Statistical Physics (1992) 

1011-1044. 

[58] D. R. M. Williams, F. MacKintosh, Driven granular media in one dimen- 
sion: Correlations and equation of state, Phys. Rev. E 54 (1996) 9-12. 

[59] Y. Zheng, H. Struchtrup, A linearization of mieussens's discrete velocity 
model for kinetic equations, Eur. J. Mech. B Fluids 26 (2007) 182-192. 



41 




t 



3 - 



Temperature profile [T(t)] for Inelastic Maxwell Molecules 
in Boltzmann Equation with Diffusion Term 



2.5 - 



H 2- 



1.5 - 



Computed 

^ Analylical 




This figure "conECorr.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/0710.5308v2 



